Aminet 33/Aminet 33 - October 1999.iso/Aminet/misc/math/TCalcStats2c.lha/TCalcStats2c/AREXX/Chi-Sq_Indep.rexx

OS/2 REXX Batch file | 1998-10-09 | 11.1 KB | 571 lines

/* Chi-Square Test for Independence rxc Contingency Table*/

options results
if ~show('P','TCALC') then do
    address command 'run turbocalc:turbocalc'
    address command 'waitforport TCALC'
    loadflag=1
end
address 'TCALC'
'DEFPUBSCREEN()'
/* Add-in Rexx Math Library needed for some routines */
signal on syntax
if ~show('l','rexxmathlib.library') then
   call addlib('rexxmathlib.library',0,-30)
if ~show('l','rexxreqtools.library') then
   call addlib('rexxreqtools.library',0,-30)
if ~show('l','rexxsupport.library') then
   call addlib('rexxsupport.library',0,-30)
 /* add to library list */
signal off syntax

/* Start Main Routine */
if loadflag=1 then 'Load()'
'ActivateWindow()'
range=rtgetstring(,"Enter Cell Range for Input","Input Request") /*,,'rt_pubscrname="TCALC"')*/
colon=pos(":",range)
if colon=0 then do
    'Message "Please select a range before executing this script"'
    'DEFPUBSCREEN("Workbench")'
    exit
end

/* Find cell references and cell, column numbers */
start_cell=substr(range,1,colon-1)
end_cell=substr(range,colon+1)
start_row=cellrow(start_cell)
end_row=cellrow(end_cell)
start_col=cellcol(start_cell)
end_col=cellcol(end_cell)
NRows=end_row-start_row+1
NCols=end_col-start_col+1

/*Labels in First Column?*/
Clabel=rtgetstring("Yes","Are there Labels in Column 1? Yes or No:","Input Request") /*,,'rt_pubscrname="TCALC"')*/
Clabel=upper(left(Clabel,1))
If Clabel="" then do
    'DEFPUBSCREEN("Workbench")'
    exit
end

/* Get cell reference for output range */
out_cell=rtgetstring(,"Enter Cell Reference for Output","Input Request") /*,,'rt_pubscrname="TCALC"')*/
if out_cell="" then do
    'DEFPUBSCREEN("Workbench")'
    exit
end
if length(out_cell)<2 | datatype(left(out_cell,1),'n')=1 then do
    'Message "Invalid cell reference"'
    'DEFPUBSCREEN("Workbench")'
    exit
end
/* Suppress Screen Redraw to Speed Things Up */
'Refresh 0'

/* Open a small output window on tcalc screen*/
fo=0
CR='0a'x
DisplayMsg="Calculating...Please Wait."||CR||"User input is disabled during calculation."||CR
if open(6Info, 'con:100/0/450/80/Progress/SCREEN TCALC', w) then do
     call writeln(6Info, DisplayMsg)
    fo=1
end
else do
    'Message "TCALC Screen not available for Progress messages"'
end
CALL DELAY(150)

/* Get cell references for top cell in each column */
top_cell.=0
'SelectCell' start_cell
do col=start_col to end_col
    'GetCursorPos'
    top_cell.col=result
    'Column 1'
end
if fo then call writech(6Info,"Progress...10 ")
col=start_col

/* Get labels for later use on output */
If Clabel="Y" then do
    'SelectCell' start_cell
    'CursorDown 1'
    'GetValue'
    do y=1 to NRows-1
        'GetValue'
        str=result
        lefttitle.y=translate(strip(str),"_"," ")
        'CursorDown 1'
    end
end
If Clabel="Y" then do
    col=start_col+1
    NCols=NCols-1
end
say "col="col
say "Clabel="Clabel
'SelectCell' top_cell.col
'GetValue'
testlabel=result
testlabel=strip(testlabel)
if datatype(testlabel,'n')=1 then do
    labelflag=0
    do x=1 to NCols
    title.x="Column_"||x
    end
end
else do
    labelflag=1
    NRows=NRows-1
    do x=1 to NCols
    'GetValue'
    str=result
    title.x=translate(strip(str),"_"," ")
    'Column 1'
    end
end
if fo then call writech(6Info,"20 ")

/* Get data from cell range */
lav=0
tot=0
count.=0
total.=0
do x=1 to NCols
    'SelectCell' top_cell.col
    if labelflag=1 then 'CursorDown 1'
    do y=1 to NRows
        'GetValue'
        valtest=result
        if datatype(valtest)='NUM' then do
            'GetValue'
            val=result
            data.x.y=val
            tot=tot+val
            total.x=tot
            count.x=1+count.x
        end
        'CursorDown 1'
    end
    col=col+1
    tot=0
    lav=0
    val=0
end
if fo then call writech(6Info,"40 ")
if count.1 ~= count.2 then do
    'Message "Equal sample sizes required - aborting!"'
    'DEFPUBSCREEN("Workbench")'
    exit
end

/* Sum of Row Totals */
Rowtot.=0
GndTot=0
Do y=1 to NRows
    Do x=1 to NCols
        Rowtot.y=(Rowtot.y)+data.x.y
    End
GndTot=GndTot+Rowtot.y
End
if fo then call writech(6Info,"60 ")

/* Calculate Expected Frequencies */

Expect.=0
Diff.=0
Chi=0
Do x=1 to NCols
    Do y=1 to NRows
        Expect.x.y=((Rowtot.y)*(total.x))/GndTot
        Diff.x.y=(((data.x.y)-(Expect.x.y))**2)/Expect.x.y
        Chi=Chi+Diff.x.y
    end
end

/* Do same with Yates Correction */
ExpectY.=0
DiffY.=0
ChiY=0
Do x=1 to NCols
    Do y=1 to NRows
        ExpectY.x.y=((Rowtot.y)*(total.x))/GndTot
        Temp=(ABS((data.x.y)-(ExpectY.x.y)))-.5
        DiffY.x.y=(Temp**2)/ExpectY.x.y
        ChiY=ChiY+DiffY.x.y
    end
end

/* Calculations for Standardised Residuals */
Res.=0
Do x=1 to NCols
    Do y=1 to NRows
        Res.x.y=((data.x.y)-(Expect.x.y))/sqrt(Expect.x.y)
    end
end

/* Measures of Association */

if NCols =2 & NRows=2 then do
    Phi=sqrt(Chi/GndTot) /* phi coefficient */
    Q1=((data.1.1)*(data.2.2))-((data.2.1)*(data.1.2))
    Q2=((data.1.1)*(data.2.2))+((data.2.1)*(data.1.2))
    YuleQ=Q1/Q2 /* Yules Q coefficient */
    Odd1=((data.1.1)/GndTot)*((data.1.2)/GndTot)
    Odd2=((data.2.1)/GndTot)*((data.2.2)/GndTot)
    OddsRatio=Odd1/Odd2 /* Odds Ratio */
end
PearC=sqrt(Chi/(Chi+GndTot)) /* Pearson contingency coeff */
Select
    when NRows <=NCols then k=NRows
    when NRows>NCols then k=NCols
otherwise
    NOP
end
PhiCr=sqrt(Chi/(GndTot*(k-1))) /* Cramers Phi Coefficient */

df=(NRows-1)*(NCols-1)

/* Calculate Probability */
P=CHIPROB(Chi,df)
PY=CHIPROB(ChiY,df)
if fo then call writech(6Info,"80 ")

/* Calculate Critical */
Chicrit1=CHICRIT(.95,df)
if fo then call writech(6Info,"90 ")
Chicrit2=CHICRIT(.99,df)

if fo then do
    call writeln(6Info,"100 ")
    call writeln(6Info,"Writing output to window...")
end

/* Output */
'SelectCell' out_cell
'ColumnWidth' 20
'Put "Chi Square Test for Independence"'
'CursorDown 2'
'Put "Expected Frequencies"'
'CursorDown 1'
'GetCursorPos'
top_cell=result
If Clabel="Y" then do
    Do y=1 to NRows
        'CursorDown 1'
        'Put' lefttitle.y
    end
'SelectCell' top_cell
'Column 1'
end
Do x=1 to NCols
    'Alignment 2'
    'Put' title.x
    'CursorDown 1'
    Do y=1 to NRows
        'Put' Expect.x.y
        'CursorDown 1'
    end
'SelectCell' top_cell
'Column 2'
end
'SelectCell' top_cell
'CursorDown 4'
'Put "Chi-Square:"'
'CursorDown 1'
'Put "d.f.:"'
'CursorDown 1'
'Put "P(CHI<=chi):"'
'CursorDown 1'
'Put "Chi-Critical (95%):"'
'CursorDown 1'
'Put "Chi-Critical (99%):"'
'CursorUp 4'
'Column 1'
'Put' format(Chi,,4)
'CursorDown 1'
'Put' df
'CursorDown 1'
'Put' format(P,,6)
'CursorDown 1'
'Put' format(ChiCrit1,,4)
'CursorDown 1'
'Put' format(ChiCrit2,,4)
'Column -1'
'CursorDown 1'
'Put "With Yates Correction:"'
'CursorDown 1'
'Put "Chi-Square:"'
'CursorDown 1'
'Put "P(CHI<=chi):"'
'CursorUp 1'
'Column 1'
'Put' format(ChiY,,4)
'CursorDown 1'
'Put' format(PY,,6)
'CursorDown 1'
'Column -1'
'Put "Measures of Association:"'
if NCols =2 & NRows=2 then do
    'CursorDown 1'
    'Column 1'
    'Put' format(Phi,,4)
    'Column -1'
    'Put "Phi Coeff.:"'
    'CursorDown 1'
    'Put "Yules Q:"'
    'Column 1'
    'Put' format(YuleQ,,4)
    'CursorDown 1'
    'Put' format(OddsRatio,,4)
    'Column -1'
    'Put "Odds Ratio:"'
end
'CursorDown 1'
'Put "Pearson Cont. Coeff.:"'
'Column 1'
'Put' format(PearC,,4)
'CursorDown 1'
'Put' format(PhiCr,,4)
'Column -1'
'Put "Cramers Phi Coeff.:"'
'Refresh 1'
'Refresh 2'
/*'Message' "Finished"*/
/*indicate the main script is finished*/
DisplayMsg="Cleaning up ...."||CR||"Exiting"
result=writeln(6Info, DisplayMsg)
if result~=0 then do
    /*Wait 3 seconds*/
    CALL DELAY(150)
    /* close window*/
    result=close(6Info)
end
'DEFPUBSCREEN("Workbench")'
exit

/* Procedures */

cellrow: procedure
do
    parse arg cell
    do charpos=2 to length(cell)
    if datatype(substr(cell,charpos,1),n) then return substr(cell,charpos)
    end
    return 0
end
Return

cellcol: procedure
do
    parse arg cell
    labels="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
    cell=upper(cell)
    len=length(cell)
    val=0
do charpos=1 to len
    if datatype(substr(cell,charpos,1),n) then
    do cell=reverse(substr(cell,1,charpos-1))
    do x=1 to length(cell)
    val=(26**(x-1))*pos(substr(cell,x,1),labels)+val
    end
    return val
    end
    end
    return 0
end
Return
/* It is important to put the exposed array at the end of the next line */
Sort: procedure expose NCols NRows data.
do x=1 to NCols
L=(xtoy(2,int(log(NRows-1)/log(2))))-1
    Do Until L<1
    L=trunc(int(L/2))
    Do J=1 to L
            Do K=J+L To NRows-1 By L
            I=K
            dumdat=data.x.I
            Do while I>L
                y=I-L
                If data.x.y ~> dumdat then Leave
                data.x.I=data.x.y
                I=I-L
            End
            data.x.I=dumdat
            End
        End
    End
End
Return

syntax:
     if arg(1)='FAIL' then do
        'Message "Library is unavailable."'
        'DEFPUBSCREEN("Workbench")'
        exit
        end
    'DEFPUBSCREEN("Workbench")'
    exit

Format:  procedure

     arg number, before, after
     CallLine = SIGL
     if ~datatype(CallLine, 'N') then CallLine = '??'

     /* Make sure we have a number as first (required) argument    */
     if ~datatype(number, 'N') then do
        if number = '' then
           fc = 17     /* Wrong number of arguments           */
        else
           fc = 47     /* Arithmetic conversion error             */
        signal FormatSyntaxError
     end
     num = number + 0
     if before = '' & after = '' then
        return num
     else do
        parse var num integer '.' fraction
        if before = '' then before = length(integer)
        if after = '' then after = length(fraction)
        if ~datatype(before, N) | ~datatype(after, N) then
           do fc = 18
           signal FormatSyntaxError
       end
        if before < length(integer) then do
           fc = 18
           signal FormatSyntaxError
        end
        if after ~= length(fraction) then do
           fraction = trunc(('.'fraction'0') + ('.'copies('0', after)'5'), after)
        if integer<1&integer>-1 then integer=integer
           else integer = integer + (fraction % 1)
           fraction = substr(fraction, 3)
        end
        if fraction >= 0 then
           return right(integer, before)'.'fraction
        else
           return right(integer, before)
     end

 FormatSyntaxError:
        if show('F', STDERR) then
           call writeln(STDERR, '+++ Error' fc 'in line' CallLine':' errortext(fc))
        else
           mess='+++ Error' fc 'in line' CallLine':' errortext(fc)
        'Message' mess
        parse source Func .
        if Func = 'FUNCTION' then do
        'DEFPUBSCREEN("Workbench")'
           exit "Err"
        end
        else do
        'DEFPUBSCREEN("Workbench")'
           exit 10
        end

CHIPROB: PROCEDURE

    PARSE ARG X2,K1
    X3=0.5*X2
    K3=0.5*K1
    X=K3+1
    LG=LOGGAMMA(X)
    S=0
    A=EXP(K3*LN(X3)-LG-X3)
    IF A~=0 THEN DO
        T=1
        S=1
        G=K3
        DO WHILE (T/S) >.000001
            G=G+1
            T=T*X3/G
            S=S+T
        END
    END
    PO=S*A
RETURN PO

CHICRIT: PROCEDURE

    PARSE ARG PO,K1
    AO=.0001
    E=.005
    E2=E+E
    PA=PO
    SELECT
    WHEN K1=1 THEN DO
        PO=.5+PA/2
        Z=NORMALPP(PO)
        X1=Z*Z
    END
    WHEN K1=2 THEN DO
        X1=1-PA
        X1=-2*LN(X1)
    END
    OTHERWISE
        Z=NORMALPP(PO)
        A1=2/(9*K1)
        W=1-A1+Z*SQRT(A1)
        XO=K1*(W**3)
        DO FOREVER
            X2=XO
            PO=CHIPROB(X2,K1)
            F1=PO-PA
            X2=XO+E
            PO=CHIPROB(X2,K1)
            F2=PO
            X2=XO-E
            PO=CHIPROB(X2,K1)
            F2=F2-PO
            F2=F2/E2
            X1=XO-F1/F2
            IF ABS(X1-XO)>AO THEN XO=X1
            ELSE LEAVE
        END
    END
    X2=X1
RETURN X2

LOGGAMMA: PROCEDURE

    ARG XA
    C1=76.18009173
    C2=-86.50532033
    C3=24.01409822
    C4=-1.231739516
    C5=.001208580
    C6=-.000005364
    C7=2.506628275
    X1=XA-1
    W=X1+5.5
    W=(X1+.5)*LN(W)-W
    S=1+C1/(X1+1)+C2/(X1+2)+C3/(X1+3)
    S=S+C4/(X1+4)+C5/(X1+5)+C6/(X1+6)
    L=W+LN(C7*S)
RETURN L

NORMALPP: PROCEDURE

    ARG P0
    A1=2.515517
    A2=.802853
    A3=.010328
    A4=1.432788
    A5=.189269
    A6=.001308
    Q=.5-ABS(P0-.5)
    SN=SIGN(-2*LN(Q))
    IF SN=-1 THEN W=-1*SQRT(ABS(-2*LN(Q))
    ELSE W=SQRT(-2*LN(Q))
    W1=A1+W*(A2+A3*W)
    W2=1+W*(A4+W*(A5+A6*W))
    ZZ=W-W1/W2
    SELECT
        WHEN (P0-.5)<0 THEN TT=-1
        WHEN (P0-.5)=0 THEN TT=0
        otherwise TT=1
    END
    ZZ=ZZ*TT
RETURN ZZ